Análisis de Datos de Área: El Analfabetismo en Nepal

Authors

Camilo Gómez

Sebastián Olarte

Descripción Del Conjunto De Datos: International Aid to Nepal (2007–2014)

El conjunto de datos utilizado corresponde al proyecto International Aid to Nepal (2007–14), desarrollado por el Center for Spatial Data Science (GeoDa Center). Su objetivo es integrar información geoespacial, socioeconómica y financiera que permita estudiar la pobreza, el desarrollo humano y la distribución de la ayuda internacional en los 75 distritos administrativos de Nepal. Este dataset combina datos espaciales y estadísticos.

Resumen General del Dataset

Característica	Detalle
Título	International Aid to Nepal (2007–14)
Observaciones	75 distritos (un polígono por distrito)
Variables	61 variables
Cobertura temporal	1991–2013 (ayuda internacional: 2007–2014)
Fuentes principales	AidData y Open Nepal
Tipo de dato espacial	Polígonos distritales (MULTIPOLYGON)
CRS	WGS84 – EPSG:4326
Archivos shapefile	`.shp`, `.shx`, `.dbf`, `.prj`, `.qpj`

Categorías y Descripción de Variables

Las 61 variables se organizan en cinco grupos temáticos principales.

A. Variables de Identificación Geográfica

Variable	Descripción	Unidad/Tipo
`id`	Identificador único del distrito	Numérica
`name_1`	Región administrativa	Texto
`name_2`	Zona administrativa	Texto
`district`	Nombre del distrito	Texto
`lon`	Longitud del centroide	Grados
`lat`	Latitud del centroide	Grados

B. Indicadores de Pobreza y Economía

Variable	Descripción	Unidad
`DepEcProv`	Privación económica	Índice
`PovIndex`	Índice de Pobreza Humana	Índice
`PCInc`	Ingreso per cápita	Rs.
`PCIncPPP`	Ingreso per cápita ajustado por PPA	USD PPP
`PCIncMP`	Ingreso per cápita a precio de mercado	Rs.

C. Indicadores Sociales, Salud y Servicios

Variable	Descripción	Unidad
`Population`	Población total	Personas
`MalKids`	% de niños malnutridos (< 5 años)	Porcentaje
`Lif40`	% que no alcanzará los 40 años	Porcentaje
`NoSafH20`	% sin acceso a agua potable segura	Porcentaje
`VotNum`	Número de votantes registrados	Personas

D. Indicadores Educativos (2012–2013)

Variable	Descripción	Unidad
`BoyG1_5`	Niños matriculados grado 1–5	Número
`GirlG1_5`	Niñas matriculadas grado 1–5	Número
`KIDS1_5`	Total niños grado 1–5	Número
`SchoolCnt`	Número de escuelas	Número
`SCHLPKID`	Escuelas por 1.000 niños	Ratio
`SCHLPPOP`	Escuelas por 1.000 habitantes	Ratio
`AD_ILLIT`	Analfabetismo adulto (2011)	Porcentaje
`AD_ILGT50`	1 si analfabetismo > 50%, 0 en otro caso	Binaria

E. Variables de Ayuda Internacional

Cada sector contiene dos variables: xxCAMT (monto comprometido) y xxDAMT (monto desembolsado).

Código	Sector	Código	Sector
AG	Agricultura	HEALT	Salud
BANK	Banca	HUM	Ayuda humanitaria
COMM	Comunicaciones	IND	Industria
CON	Conflicto	MUL	Multisectorial
EDU	Educación	SOC	Infraestructura social
ENGY	Energía	TOUR	Turismo
ENV	Medio ambiente	TRAN	Transporte y almacenamiento
FOR	Silvicultura	WAT	Agua y saneamiento
GOV	Gobierno y sociedad civil	TOT	Total general

Lectura de Datos

Se realiza la lectura de los datos de Nepal.

Code

nepal <- st_read("nepal/nepal.shp", quiet = TRUE) # Se leen los datos de Nepal

nepal.st <- st_geometry(nepal) # Polígonos

#st_crs(nepal) # CRS

# Transformación a UTM

nepal.utm <- st_transform(nepal, crs = 32645)

La Figure 1 muestra el mapa coroplético de la tasa de analfabetismo adulto en Nepal a nivel distrital. El mapa evidencia una marcada heterogeneidad espacial, con valores elevados concentrados principalmente en las regiones occidental y noroccidental del país, así como en parte de la franja sur, mientras que los niveles más bajos se observan en los distritos del centro y este.

La distribución espacial de los valores sugiere la presencia de agrupamientos de distritos con tasas similares de analfabetismo, indicando que la variable no se distribuye de manera aleatoria en el territorio. Este patrón visual constituye una evidencia descriptiva inicial que motiva la aplicación de métodos de análisis espacial para evaluar formalmente la dependencia espacial.

tmap_mode("view")

tm_shape(nepal.utm) +
  tm_polygons("ad_illit",
              fill.scale = tm_scale_intervals(n = 5, 
                                              style = "quantile", 
                                              values = "brewer.reds"),
              fill.legend = tm_legend(title = "Tasa de\nAnalfabetismo")) +
  tm_title("") +
  tm_layout(legend.position = c("left", "bottom"), frame = FALSE) +
  tm_scalebar(position = c("right", "bottom")) +
  tm_compass(position = c("right", "top"))

Figure 1: Mapa de la tasa de analfabetismo en Nepal por distrito

De esta forma, es posible ver que hay agrupaciones de distritos cuya tasa de analfabetismo es alta y otros para la cual la tasa es baja, indicando presencia de correlación espacial. Hacia el noroccidente y suroriente de Nepal se encuentran los distritos con una mayor tasa de analfabetismo.

Matriz de Vecindades

El primer paso consiste en definir los centroides de las áreas, que en este caso corresponden a los distritos de Nepal. Para ello, se emplean los centroides geométricos de cada distrito.

Code

coords <- st_coordinates(st_centroid(nepal.utm))

Pesos Espaciales

De esta forma, se procede con la definición de los vecinos para evaluar cuál es la matriz que mejor describe la relación espacial presente en los datos de analfabetismo de Nepal. En el caso de los vecinos físicos se obtienen las mismas conexiones para los métodos queen y rook.

Code

## Vecinos Físicos

nepal.nb <- poly2nb(nepal.utm, queen = TRUE)
nepal.lw_w <- nb2listw(nepal.nb, style = "W", zero.policy = TRUE)
nepal.lw_b <- nb2listw(nepal.nb, style = "B", zero.policy = TRUE)

## Vecinos basados en grafos

tri.nb <- tri2nb(coords)
tri.nb_w <- nb2listw(tri.nb, style = "W", zero.policy = TRUE)
tri.nb_b <- nb2listw(tri.nb, style = "B", zero.policy = TRUE)


soi.nb <- graph2nb(soi.graph(tri.nb, coords))
soi.nb_w <- nb2listw(soi.nb, style = "W", zero.policy = TRUE)
soi.nb_b <- nb2listw(soi.nb, style = "B", zero.policy = TRUE)


relative.nb <- graph2nb(relativeneigh(coords))
relative.nb_w <- nb2listw(relative.nb, style = "W", zero.policy = TRUE)
relative.nb_b <- nb2listw(relative.nb, style = "B", zero.policy = TRUE)


gabriel.nb <- graph2nb(gabrielneigh(coords), sym=TRUE)
gabriel.nb_w <- nb2listw(gabriel.nb, style = "W", zero.policy = TRUE)
gabriel.nb_b <- nb2listw(gabriel.nb, style = "B", zero.policy = TRUE)


## Vecinos basados en distancia

knn.nb_1 <- knn2nb(knearneigh(coords, k = 1))
knn.nb_2 <- knn2nb(knearneigh(coords, k = 2))
knn.nb_3 <- knn2nb(knearneigh(coords, k = 3))
knn.nb_4 <- knn2nb(knearneigh(coords, k = 4))


knn.nb_1_w <- nb2listw(knn.nb_1, style = "W", zero.policy = TRUE)
knn.nb_1_b <- nb2listw(knn.nb_1, style = "B", zero.policy = TRUE)


knn.nb_2_w <- nb2listw(knn.nb_2, style = "W", zero.policy = TRUE)
knn.nb_2_b <- nb2listw(knn.nb_2, style = "B", zero.policy = TRUE)


knn.nb_3_w <- nb2listw(knn.nb_3, style = "W", zero.policy = TRUE)
knn.nb_3_b <- nb2listw(knn.nb_3, style = "B", zero.policy = TRUE)


knn.nb_4_w <- nb2listw(knn.nb_4, style = "W", zero.policy = TRUE)
knn.nb_4_b <- nb2listw(knn.nb_4, style = "B", zero.policy = TRUE)

Luego de haber definido las matrices de conectividad, se realiza la elección de la mejor con ayuda del test I de Morán: se realiza el test y se escoge la que de un menor \(p\)-valor. Primero se almacenan todas las matrices en una lista y se define una función para realizar el procedimiento.

El sistema de hipótesis que se quiere juzgar es el siguiente, donde \(I\) representa el índice de Morán (Moraga (2023)):

\[ \begin{cases} H_0: I= \mathrm{E}[I] & \text{(no hay autocorrelación espacial)}\\ \text{vs} \\ H_1: I\not=\mathrm{E}[I] & \text{(hay autocorrelación espacial)}\\ \end{cases} \]

mat <- list(nepal.lw_w, nepal.lw_b,
            tri.nb_w, tri.nb_b,
            soi.nb_w, soi.nb_b,
            relative.nb_w, relative.nb_b,
            gabriel.nb_w, gabriel.nb_b,
            knn.nb_1_w, knn.nb_1_b,
            knn.nb_2_w, knn.nb_2_b,
            knn.nb_3_w, knn.nb_3_b,
            knn.nb_4_w, knn.nb_4_b)


names(mat) <- c(paste0("phys_", c("W","B")),
                paste0("tri_",  c("W","B")),
                paste0("soi_",  c("W","B")),
                paste0("rel_",  c("W","B")),
                paste0("gab_",  c("W","B")),
                paste0("knn1_", c("W","B")),
                paste0("knn2_", c("W","B")),
                paste0("knn3_", c("W","B")),
                paste0("knn4_", c("W","B")))

matrix_eval <- function(var, weights) {
  
  aux <- numeric(length(weights)) 
  
  for (i in seq_along(weights)) {
    aux[i] <- moran.mc(var,
                       weights[[i]], nsim = 9999,
                       alternative = "two.sided")$p.value
  }
  
  best_idx <- which.min(aux)
  
  list(
    best_index  = best_idx,
    best_weight = weights[[best_idx]],
    best_name = names(weights)[which.min(aux)],
    best_moran  = moran.mc(var, weights[[best_idx]], nsim = 9999,
                           alternative = "two.sided"),
    all_pvalues = aux
  )
}

Se ve cuál fue la matriz que mejor describe la correlación espacial de los datos.

Code

best.w <- matrix_eval(nepal.utm$ad_illit, mat)
best.weights <- best.w$best_weight

best.weights

Characteristics of weights list object:
Neighbour list object:
Number of regions: 75 
Number of nonzero links: 366 
Percentage nonzero weights: 6.506667 
Average number of links: 4.88 

Weights style: W 
Weights constants summary:
   n   nn S0       S1       S2
W 75 5625 75 32.85108 309.2524

Como se puede ver, esta corresponde a la obtenida considerando los cuatro vecinos más cercanos con el estilo W, es decir, con las filas estandarizadas. El mapa con esta matriz de conectividad se puede ver en la Figure 2.

nb_lines <- nb2lines(best.weights$neighbours, coords = coords, as_sf = TRUE)

st_crs(nb_lines) <- st_crs(nepal.utm)

ggplot() +
  geom_sf(data = nepal.utm, fill = "grey95", color = "grey80", size = 0.2) +
  # Conexiones
  geom_sf(data = nb_lines, color = "#2c3e50", size = 0.3, alpha = 0.6) +
  
  # Centroides
  geom_point(data = coords, aes(x = X, y = Y), 
             color = "#e74c3c", size = 1, alpha = 0.8) +
  
  # Estética
  theme_void() +
  labs(title = "Estructura de Conectividad Espacial",
       subtitle = paste0("Matriz de Pesos: ", best.w$best_name))

Figure 2: Mapa de Nepal con la Mejor Estructura de Conectividad Espacial

Se realizó el Test I de Moran mediante Monte Carlo con \(n=9999\) permutaciones. La distribución empírica del índice de Moran puede observarse en la Figure 3.

Figure 3: Histograma de los valores simulados del índice de Moran. La línea cortada representa del valor del índice de Moran observado.

De esta forma, se evidencia que al 5% de significancia se rechaza \(H_0: I = \mathrm{E}[I]\). Es decir, hay suficiente evidencia para concluir que hay presencia de correlación espacial con respecto a la tasa de adultos analfabetas en Nepal.

En la Figure 4 se muestra el diagrama de Moran de la tasa de adultos analfabetas en Nepal.

Figure 4: Dispersograma de Moran para la tasa de adultos analfabetas en Nepal

Índice de Moran (\(I = 0.595\)):
El valor del estadístico es positivo y sustancialmente mayor que cero, lo que indica la presencia de una autocorrelación espacial global positiva fuerte*. Esto sugiere que los distritos geográficamente cercanos tienden a presentar tasas de analfabetismo similares: áreas con niveles altos tienden a agruparse espacialmente, al igual que áreas con niveles bajos.
Dado que el \(p\)-valor es prácticamente cero (\(p < 0.001\)), se rechaza la hipótesis nula de ausencia de autocorrelación espacial. En consecuencia, el patrón observado no es producto del azar, sino que refleja una estructura espacial estadísticamente significativa.
1. Cuadrantes

El diagrama se divide en cuatro cuadrantes definidos por las medias de la tasa de analfabetismo (eje horizontal) y su rezago espacial (eje vertical). Cada cuadrante representa un tipo particular de asociación espacial:

Cuadrante	Patrón	Descripción	Ejemplos en la Figura	Implicación Geográfica
Arriba–Derecha	Alto–Alto (H–H)	Distritos con tasas de analfabetismo altas rodeados por vecinos con valores igualmente altos.	Casos en el extremo superior derecho.	Hotspots de analfabetismo.
Abajo–Izquierda	Bajo–Bajo (L–L)	Distritos con tasas de analfabetismo bajas rodeados por vecinos con valores bajos.	Lalitpur y puntos cercanos al origen.	Coldspots de analfabetismo.
Arriba–Izquierda	Bajo–Alto (L–H)	Distritos con tasas bajas rodeados por vecinos con tasas altas.	Surkhet.	Outliers espaciales (rompen el patrón regional).
Abajo–Derecha	Alto–Bajo (H–L)	Distritos con tasas altas rodeados por vecinos con tasas bajas.	Kapilbastu.	Outliers espaciales con comportamiento diferenciado.

La pendiente positiva de la recta de regresión —equivalente al valor del Índice de Moran— y la concentración de observaciones en los cuadrantes Alto–Alto y Bajo–Bajo confirman que la tasa de analfabetismo adulto es un fenómeno espacialmente dependiente en Nepal. Este resultado sugiere que factores subyacentes al analfabetismo (por ejemplo, inversión en educación, pobreza estructural o condiciones socioculturales) no se distribuyen de manera aleatoria, sino que operan a escala regional y trascienden los límites administrativos distritales. En consecuencia, el diseño de políticas públicas orientadas a la reducción del analfabetismo debería adoptar un enfoque regional o interdistrital, concentrando esfuerzos en los hotspots identificados y analizando los outliers espaciales para comprender mecanismos locales y replicar experiencias exitosas.

Prueba de Autocorrelación Local

Con el fin de identificar clusters y evaluar la similitud entre áreas y sus vecinos de forma local, se realiza el test I de Moran de forma locaL. En la Figure 5 se muestran dos mapas en los que se observa autocorrelación espacial significativa en los distritos de las zonas noroccidental, central y suroriental de Nepal. En estas zonas predominan relaciones tipo alto–alto y bajo–bajo, lo que indica que los distritos con altas tasas de analfabetismo en adultos tienden a estar rodeados por distritos con tasas igualmente altas, y de igual forma para las tasas bajas.

lmoran <- localmoran(nepal.utm$ad_illit, 
                     best.weights, 
                     alternative = "two.sided")
nepal.utm$lmp <- lmoran[, "Pr(z != E(Ii))"]

Figure 5: Análisis de Autocorrelación Espacial Local para la tasa de analfabetismo.

Modelos de Regresión

Se ajustan distintos modelos de regresión, teniendo en cuenta las variables descritas al inicio del documento.

data <- as.data.frame(nepal.utm)

data <- data %>% 
  mutate(ag = AGDAMT / AGCAMT,
         engy = log((ENGYCAMT + 1 )/ (ENGYDAMT + 1)),
         gov = (GOVDAMT + 1) / (GOVCAMT + 1),
         healt = HEALTDAMT / population,
         mult = MULTDAMT / MULTCAMT,
         soc = (SOCDAMT + 0.001) / (SOCCAMT + 0.001),
         tot = TOTDAMT / population,
         log.schlppop = log(schlppop))

Modelo de Regresión MCO

Primero, se muestra el modelo de regresión que tuvo un mejor ajuste con respecto al \(R^2\) (0.4917) y la significancia de los coeficientes. Este modelo tiene la forma \[ \mathbf{y} = \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\varepsilon}, \qquad \boldsymbol{\varepsilon} \sim \mathrm{MVN}(\mathbf{0}, \sigma^2 \mathbf{I}) \]

Code

form <- as.formula("ad_illit ~ pcinc + boyg1_5 + lif40 + soc + tot")

fit.OLS <- lm(form,
              data = data)

	Estimate	Std. Error	t value	Pr(>\|t\|)
(Intercept)	57.9031	11.9913	4.8288	0.0000
pcinc	-0.0231	0.0043	-5.4282	0.0000
boyg1_5	0.0002	0.0001	3.5416	0.0007
lif40	1.5015	0.5215	2.8790	0.0053
soc	-25.6299	12.2164	-2.0980	0.0396
tot	0.0096	0.0049	1.9531	0.0549

Además, se realiza el Test I de Moran con los residuales del modelo, rechazando la ausencia de autocorrelación espacial al 5% de significancia. Es decir, los residuales presentan una dependencia espacial que no se pudo explicar con el modelo.

lmmorantest <- lm.morantest(fit.OLS, listw = best.weights)
t.lmmorantest <- tidy(lmmorantest)

kable(t.lmmorantest[,c("statistic", "p.value", "method", "alternative")])

statistic	p.value	method	alternative
6.216516	0	Global Moran I for regression residuals	greater

De esta forma, se hace necesario ajustar modelos de regresión espacial.

SAR - Spatial Autoregressive Model (Spatial Lag)

El modelo SAR incorpora la interdependencia espacial de la variable espacial con el \(\rho \in [-1,1]\). Este modelo tiene la estructura: \[ \mathbf{y} = \rho \mathbf{W} \mathbf{y} + \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\varepsilon}, \qquad \boldsymbol{\varepsilon} \sim \mathrm{MVN}(\mathbf{0}, \sigma^2 \mathbf{I}) \]

Code

fit.SAR <- lagsarlm(form,
                    data = data, listw = best.weights)

SEM - Spatial Error Model

Para el modelo SEM se tiene la estructura \[ \begin{aligned} \mathbf{y} =& \mathbf{X} \boldsymbol{\beta} + \mathbf{u}, \\ \mathbf{u} =& \lambda \mathbf{W} \mathbf{u} + \boldsymbol{\varepsilon}, \qquad \boldsymbol{\varepsilon} \sim \mathrm{MVN}(\mathbf{0}, \sigma^2 \mathbf{I}) \end{aligned} \] y este asume que hay factores no observados que generan la dependencia espacial. Es decir, la dependencia se incorpora en los errores.

Code

fit.SEM <- errorsarlm(form,
                      data = data, listw = best.weights)

SDM - Spatial Durbin Model

El modelo de Durbin integra tanto interdependencia espacial de la respuesta como efectos espaciales de las covariables. \[ \mathbf{y} = \rho \mathbf{W} \mathbf{y} + \mathbf{X} \boldsymbol{\beta} + \gamma \mathbf{W} \mathbf{X} + \boldsymbol{\varepsilon} , \qquad \boldsymbol{\varepsilon} \sim \mathrm{MVN}(\mathbf{0}, \sigma^2 \mathbf{I}) \]

Code

fit.Durbin <- lagsarlm(form,
                       data = data, listw = best.weights, type = "mixed")

SLX - Spatially Lagged X Model

En el modelo SLX se tienen en cuenta las relaciones entre las covariables de los vecinos. \[ \mathbf{y} = \mathbf{X} \boldsymbol{\beta} + \gamma \mathbf{W} \mathbf{X} + \boldsymbol{\varepsilon} , \qquad \boldsymbol{\varepsilon} \sim \mathrm{MVN}(\mathbf{0}, \sigma^2 \mathbf{I}) \]

Code

fit.D_lagX <- lmSLX(form,
                    data = data, listw = best.weights)

SDEM - Spatial Durbin Error Model

El modelo SDEM combina los modelos SEM y SLX. Este tiene la siguiente estructura \[ \begin{aligned} \mathbf{y} =& \mathbf{X} \boldsymbol{\beta} + \lambda \mathbf{W} \mathbf{X} + \mathbf{u}\\ \mathbf{u} =& \gamma\mathbf{W} \mathbf{u} + \boldsymbol{\varepsilon}, \qquad \boldsymbol{\varepsilon} \sim \mathrm{MVN}(\mathbf{0}, \sigma^2 \mathbf{I}) \end{aligned} \]

Code

fit.SDEM <- errorsarlm(form,
                       data = data, listw = best.weights, etype = "emixed")

Manski Model

Es un modelo complejo de la forma \[ \begin{aligned} \mathbf{y} =& \rho \mathbf{W} \mathbf{y} + \mathbf{X} \boldsymbol{\beta} + \mathbf{W} \mathbf{X} \theta + \mathbf{u}, \\ \mathbf{u} =& \lambda \mathbf{W} \mathbf{u} + \boldsymbol{\varepsilon}, \qquad \boldsymbol{\varepsilon} \sim \mathrm{MVN}(\mathbf{0}, \sigma^2 \mathbf{I}) \end{aligned} \]

Code

fit.Manski <- sacsarlm(form,
                       data = data, listw = best.weights, type="sacmixed")

SAC (SARAR)- Combined Spatial Autocorrelation Model

Por último, se tiene una combinación de modelos SAR y SEM. Es decir, se asume interdependencia espacial de la variable dependiente y factores no observados. \[ \begin{aligned} \mathbf{y} =& \rho \mathbf{W} \mathbf{y} + \mathbf{X} \boldsymbol{\beta} + \mathbf{u} \\ \mathbf{u} =& \lambda \mathbf{W} \mathbf{u}+ \boldsymbol{\varepsilon}, \qquad \boldsymbol{\varepsilon} \sim \mathrm{MVN}(\mathbf{0}, \sigma^2 \mathbf{I}) \end{aligned} \]

Code

fit.SARAR <- sacsarlm(form,
                      data = data, listw = best.weights, type="sac")

Selección del Modelo

Para realizar la selección del mejor modelo, se recurre a comparar los AIC.

Table 1: Comparación del AIC de los modelos ajustados

Modelo	gl	AIC
fit.SARAR	9	499.7997
fit.SDEM	13	501.6684
fit.Durbin	13	501.7228
fit.SEM	8	503.7218
fit.SAR	8	504.1058
fit.D_lagX	12	531.7039
fit.OLS	7	539.9456

En la Table 1 se puede observar que el modelo con un menor AIC corresponde al SARAR.

Interpretación

Debido a que el modelo SAC/SARAR fue el de menor AIC, se muestra su resúmen. Primero, se puede ver que no todas las variables son significativas, a diferencia del modelo de regresión ajustado por mínimos cuadrados ordinarios.

coef(summary(fit.SARAR)) %>% 
  kable(format = "html", digits = 4, align = 'c', escape = FALSE)

	Estimate	Std. Error	z value	Pr(>\|z\|)
(Intercept)	32.4605	11.3319	2.8645	0.0042
pcinc	-0.0183	0.0032	-5.6597	0.0000
boyg1_5	0.0001	0.0000	1.6876	0.0915
lif40	0.9806	0.4439	2.2090	0.0272
soc	-10.1646	8.2507	-1.2320	0.2180
tot	0.0059	0.0036	1.6328	0.1025

Por otro lado, se presentan las estimaciones los parámetros \(\rho\) y \(\lambda\).

params.sac <- data.frame(Parámetro = c("$\\rho$", "$\\lambda$"),
                         Estimación = c(fit.SARAR$rho, fit.SARAR$lambda),
                         SE = c(fit.SARAR$rho.se, fit.SARAR$lambda.se))

params.sac %>% 
  kable(format = "pandoc", digits = 4, align = 'c', row.names = FALSE)

Parámetro	Estimación	SE
\(\rho\)	0.4682	0.1891
\(\lambda\)	0.5683	0.2027

De esta forma, con \(\hat{\rho}=0.468\) se confirma que hay dependencia espacial con respecto a la tasa de adultos analfabetas en Nepal. Más concretamente, un distrito puede afectar ``positivamente’’ al valor de sus vecinos, justo como se veía en la Figure 5. Por otro lado, \(\hat{\lambda} = 0.568\) indica que hay factores no incluídos en el modelo que presentan dependencia espacial. En este caso, ambas estimaciones son significativas al 5%.

Por lo tanto, se podría ver que a medida que los ingresos per cápita aumentan en un distrito, su tasa de analfabetismo disminuye. Por otro lado, si en un distrito el porcentaje de personas que no sobreviven hasta los 40 años aumenta, se espera que la tasa de analfabetismo aumente también.

Validación

Correlación Espacial

La primera evaluación que se realiza corresponde a ver si los residuales del mejor modelo ajustado siguen presentando autocorrelación espacial o no. Para esto, se realiza el Test I de Moran, concluyendo que, a una significancia del 5%, no se rechaza la hipótesis de ausencia de autocorrelación espacial de los residuales.

data$residuals <- resid(fit.SARAR)

moran.sac <- moran.test(data$residuals, listw = best.weights, alternative = "two.sided")
kable(tidy(moran.sac)[,c("statistic", "p.value", "method", "alternative")])

statistic	p.value	method	alternative
0.2946557	0.7682569	Moran I test under randomisation	two.sided

Homoscedasticidad

El siguiente supuesto que se prueba es el de varianza homogénea. Primero, en la Figure 6 se puede ver que los residuales no tienen ningún comportamiento que indique heterogeneidad de la varianza.

Figure 6: Residuales vs Valores Ajustados para el modelo SARAR

Además, se usa el test de Breusch-Pagan para juzgar \(H_0: \text{la varianza es homogénea}\).

bp.sac <- bptest.Sarlm(fit.SARAR)
kable(tidy(bp.sac)[,c("statistic", "p.value", "method")])

statistic	p.value	method
3.93796	0.5583824	studentized Breusch-Pagan test

Así, al 5% de significancia no se rechaza la hipótesis nula de homoscedasticidad.

Normalidad

En la Figure 7 se puede ver que los residuales tienden a estar cerca de la línea teórica, indicando un ajuste bueno.

Figure 7: QQ-plot de los residuales del modelo SARAR

Para corroborar esto, se realiza el test de Anderson–Darling, concluyendo que no se rechaza la hipótesis de normalidad de los residuales al 5% de significancia.

ad.sarar <- ad.test(data$residuals)
kable(tidy(ad.sarar))

statistic	p.value	method
0.6159227	0.105192	Anderson-Darling normality test

Estimación

A continuación, en la Figure 8 se presentan las estimaciones de la tasa de analfabetismo en adultos en Nepal.

nepal.utm$fitted <- fitted(fit.SARAR)
nepal.utm$residuals <- resid(fit.SARAR)

Figure 8: Mapa de predicciones de la tasa de analfabetismo en Nepal por distrito

Otros Modelos

Se decide comparar los resultados con otros modelos que no incorporan parámetros de asociación espacial como se vió en la sección anterior.

Moran Eigenvalues

Para realizar el uso de los vectores propios de la matriz \(\mathbf{M}\mathbf{W}\mathbf{M}\), con \(\mathbf{M} = \mathbf{I} - \mathbf{X}(\mathbf{X}^t\mathbf{X})^{-1}\mathbf{X}^t\), se usa la función SpatialFiltering() que toma el subconjunto de vectores propios que reducen la correlación espacial de los residuales en el modelo de regresión MCO.

Code

SF <- SpatialFiltering(form, data = data, nb = best.weights$neighbours, 
                       style = "W", verbose = FALSE, ExactEV = TRUE)

fit.SF <- lm(ad_illit ~ pcinc + boyg1_5 + lif40 + soc + tot + fitted(SF),
             data = data)

Por lo tanto, también se presenta el resumen de este modelo, que además tiene un \(R^2=0.8076\).

	Estimate	Std. Error	t value	Pr(>\|t\|)
(Intercept)	57.9031	7.1238	8.1281	0.0000
pcinc	-0.0231	0.0025	-9.1371	0.0000
boyg1_5	0.0002	0.0000	5.9614	0.0000
lif40	1.5015	0.3098	4.8462	0.0000
soc	-25.6299	7.2576	-3.5315	0.0008
tot	0.0096	0.0029	3.2876	0.0017
fitted(SF)vec3	-38.4032	4.9944	-7.6892	0.0000
fitted(SF)vec11	-21.0116	4.9944	-4.2070	0.0001
fitted(SF)vec15	21.8552	4.9944	4.3759	0.0001
fitted(SF)vec1	12.0327	4.9944	2.4092	0.0193
fitted(SF)vec4	10.6444	4.9944	2.1313	0.0375
fitted(SF)vec2	-9.9101	4.9944	-1.9842	0.0521
fitted(SF)vec10	-10.3387	4.9944	-2.0701	0.0431
fitted(SF)vec12	-9.5249	4.9944	-1.9071	0.0616
fitted(SF)vec21	-13.9456	4.9944	-2.7922	0.0071
fitted(SF)vec19	11.4086	4.9944	2.2843	0.0262
fitted(SF)vec8	8.1828	4.9944	1.6384	0.1069
fitted(SF)vec20	-9.6081	4.9944	-1.9238	0.0595
fitted(SF)vec9	-6.9303	4.9944	-1.3876	0.1708

SF.moran <- lm.morantest(fit.SF, best.weights, alternative = "two.sided")
t.SFmoran <- tidy(SF.moran)

kable(t.SFmoran[,c("statistic", "p.value", "method", "alternative")])

statistic	p.value	method	alternative
-0.0167557	0.9866315	Global Moran I for regression residuals	two.sided

Code

nepal.utm$residuals.SF <- resid(fit.SF)

SF.bptest <- bptest(fit.SF)

scatter.smooth(resid(fit.SF) ~ fitted(fit.SF),
               xlab = "Valores Ajustados", ylab = "Residuales")
legend("bottomright", 
       legend=c(paste0("Breusch-Pagan:", round(SF.bptest$statistic, 4)), 
                paste0(expression(p), "-valor: ", round(SF.bptest$p.value, 4))),
       cex=1,
       bg='salmon')

Code

SF.normal <- ad.test(resid(fit.SF))

qqnorm(resid(fit.SF),
       xlab = "Cuantiles Teóricos",
       ylab = "Cuantiles Muestrales", 
       main = "QQ-plot de los residuales")
qqline(resid(fit.SF))

legend("bottomright", 
       legend=c(paste0("Anderson-Darling:", round(SF.normal$statistic, 4)), 
                paste0(expression(p), "-valor: ", round(SF.normal$p.value, 4))),
       cex=1,
       bg='salmon')

De esta forma, es posible ver que, por ejemplo, los residuales ya no presentan correlación espacial y el supuesto de normalidad se cumple. Lo único que no se cumple es la homoscedasticidad.

GAM

Otra opción es usar los modelos GAM. Una ventaja de ellos, es que permite modelar covariables que no tengan una relación lineal con la variable respuesta. La forma de incluir la correlación espacial es mediante las coordenadas: s(X, Y), como se muestra en Bivand, Pebesma, and Gómez-Rubio (2008).

Code

library(mgcv)

data <- cbind(data, coords)

fit.gam <- gam(ad_illit ~ s(pcinc, bs = "cr") + boyg1_5
                          + s(lif40, bs = "cr") + s(soc, bs = "cr")
                          + s(tot, bs = "cr") + s(X, Y), data = data)

La única forma de interpretar los suavizamientos es con los efectos parciales de las covariables sobre la variable respueta. De la misma forma, se presenta su diagnóstivco.

nepal.utm$residuals.gam <- resid(fit.gam)

gam.moran <-  moran.test(nepal.utm$residuals.gam, best.weights, alternative =  "two.sided")

tgam.moran <- tidy(gam.moran)

kable(tgam.moran[,c("statistic", "p.value", "method", "alternative")])

statistic	p.value	method	alternative
-2.492309	0.0126916	Moran I test under randomisation	two.sided

Code

gam.bptest <- bptest(fit.gam)

scatter.smooth(resid(fit.gam) ~ fitted(fit.gam),
               xlab = "Valores Ajustados", ylab = "Residuales")
legend("bottomright", 
       legend=c(paste0("Breusch-Pagan:", round(gam.bptest$statistic, 4)), 
                paste0(expression(p), "-valor: ", round(gam.bptest$p.value, 4))),
       cex=1,
       bg='salmon')

Code

gam.normal <- ad.test(resid(fit.gam))

qqnorm(resid(fit.gam),
       xlab = "Cuantiles Teóricos",
       ylab = "Cuantiles Muestrales", 
       main = "QQ-plot de los residuales")
qqline(resid(fit.gam))

legend("bottomright", 
       legend=c(paste0("Anderson-Darling:", round(gam.normal$statistic, 4)), 
                paste0(expression(p), "-valor: ", round(gam.normal$p.value, 4))),
       cex=1,
       bg='salmon')

El supuesto de no correlación espacial de los residuales se rechaza con una significancia del 5%. Con la introducción de las coordenadas no se alcanza a elimianar la dependencia espacial de los residuales.

Comparación Predicciones

Se presentan las predicciones para los tres modelos ajustados en la Figure 10.

nepal.utm$fitted.SF <- fitted(fit.SF)
nepal.utm$fitted.gam <- fitted(fit.gam)

Figure 10: Comparación de las de predicciones de la tasa de analfabetismo en Nepal por distrito

References

Bivand, Roger S., Edzer J. Pebesma, and Virgilio Gómez-Rubio. 2008. Applied Spatial Data Analysis with r. Use r! Springer. https://doi.org/10.1007/978-0-387-78171-6.

Bohórquez, Martha. 2024. “Estadística Espacial y Espacio-Temporal Para Campos Aleatorios Escalares y Funcionales: Notas de Clase.” Universidad Nacional de Colombia.

Center for Spatial Data Science (CSDS). 2017. “International Aid to Nepal (2007–14).” https://geodacenter.github.io/data-and-lab/nepal/.

Moraga, Paula. 2023. Spatial Statistics for Data Science: Theory and Practice with r. Chapman & Hall/CRC Data Science Series. Chapman & Hall/CRC.

Rüttenauer, Tobias. 2024. “Spatial Data Analysis.” https://arxiv.org/abs/2402.09895.